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We calculate the entanglement entropy using a SU(3) quenched lattice gauge simulation. We find 
that the entanglement entropy scales as 1 /I 2 at small I as in the conformal field theory. Here / is 
the size of the system, whose degrees of freedom is left after the other part are traced out. The 
derivative of the entanglement entropy with respect to I hits zero at about I* = 0.6 ~ 0.7 [fm] and 
vanishes above the length. It may imply that the Yang-Mills theory has the mass gap of the order 
of 1//*. Within our statistical errors, no discontinuous change can be seen in the entanglement 
entropy. 

We discuss also a subtle point appearing in gauge systems when we divide a system with cuts. 
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1. Introduction 

Entanglement properties of quantum systems have been received much attention in quantum 
information theory and condensed matter physics. A simple system composed of two spin- 1/2 
system in spin singlet state is a typical example of the entangled state. Entanglement entropy is 
one of quantities measuring quantum entanglement. It can be defined in any quantum systems, 
including quantum mechanical systems and quantum field theories. The entanglement entropy 
(also called geometric entropy) between two regions, a subregion A of size / and its complement 
B, measures how much two regions are quantumly correlated and it is expected to be very useful to 
investigate phase structures of quantum system. 

Quantum entanglement of ground states has been widely studied in condensed matter physics 
(for a review, see [[]]]). In the Ising chain model, for instance, the entanglement entropy at the 
critical point diverges while it saturates in the non-critical regime. The entanglement entropy can 
serve as an order parameter of quantum phase transitions. 

As stated above, entanglement entropy can be defined in quantum field theories. The pure 
Yang-Mills theory is particularly interesting since it is a confining theory and is expected to have 
a mass gap. Recently, gauge/gravity duality has been extensively studied and it provides a method 
to study non-perturbative infrared dynamics of confining gauge theories. The calculation of the 
entanglement entropy using holographic approach has been proposed by Ryu and Takayanagi [Q] 
(for a review on the holographic calculation, see [Q] ), and generalized by Klebanov et al.|Q], as 
the minimal surface 7 bending down to the bulk space, 



Here is the 10 dimensional Newton constant, G- nc | the induced string frame metric on the 
surface 7, and the dilaton field. Although gravitational background dual to pure Yang-Mills 
theory has not been discovered, some approaches have been proposed including an effective model, 
so-called AdS/QCD. The numerical simulations of the entanglement entropy will give valuable 
numerical support to this kind of holographic approach. 

In the holographic approach, the calculation of the entanglement entropy in the gauge theory 
side is reduced to the calculation of geodesies in the gravity side. The boundary of geodesies 
coincides the boundary of partitioned subsystems. This is quite similar to the calculations of the 
Wilson loop in the holographic approach. On the gravity side, Wilson loops are obtained by the 
action of the string world sheet whose boundary is the Wilson loop. The entanglement entropy 
has been studied for various confining backgrounds [Q, f50. For some of confining backgrounds, 
two surfaces, called connected and disconnected surfaces, compete and the former dominates at 
small /. At some critical length I*, disconnected surface dominates and the entanglement entropy 
becomes /-independent at large /. The expected behavior of the entanglement entropy for the AdS 
bubble solution is schematically depicted in Fig. [|. The connected surface is the order of N% in 
large N expansion while disconnected surfaces are 0(1). This indicates that the effective degrees 
of freedom at small / are gluonic degrees of freedom. By contrast, those at large I are glueballs, 
color singlet objects. Therefore, it may be natural to expect that the critical length I* plays the role 
of (the inverse of) the critical temperature T c of the confinement/deconfinement phase transition. 
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Figure 1: Schematic picture of the entanglement entropy predicted by the holographic approach for the 
AdS bubble solution. At small I, 3Sa/31 behaves as 1 /I 3 as conformal field theories in (3 + 1) -dimensional 
spacetime. By contrast, it vanishes at large I where disconnected surfaces dominate. 

The entanglement entropy in SU(2) lattice gauge theory has been studied by Velytsky ^ and 
Bividovich and Polikarpov n/n. In Ref.[|]], SU(N) lattice gauge theories are studied in Migdal- 
Kadanoff approximation, and in Ref.[|7]], SU(2) lattice gauge theory is numerically investigated, 
and there is an indication that the derivative of the entanglement entropy shows a discontinuous 
change at some critical length scale I* and it vanishes. 

In this paper, we investigate the entanglement entropy in SU(3) pure Yang-Mills theory using 
lattice Monte Carlo simulations. Instead of directly calculating the entropy, we adopt numerical 
technique to evaluate the entanglement entropy, which has also been used in [0] (originally pro- 
posed in [|8[ |9|] in order to calculate the pressure in the deconfined phase). 



2. Definition and properties of the entanglement entropy 




Figure 2: The complementary regions A 
and B separated by an imaginary boundary 
at x = I. y and z axes are perpendicular to 
the plane. Separation is purely an imag- 
inary process and nothing has to be done 
on the physical state. The entanglement 
entropy measures quantum correlation be- 
tween two regions A and B. 



The entanglement entropy of a pure state is 
defined as follows. We divide the total system into sub- 
region A and its complement B. See Fig. ||. Let / be the 
size of the system A in the x direction. The density ma- 
trix of the system is p = |*F) . Since we consider the 
pure state, the von Neumann entropy of the system is 
clearly zero. The reduced density matrix obtained by 
tracing out the degrees of freedom in the region B, 



p A =Tv B p=Tr B \^)^\ 



(2.1) 



describes the density matrix for an observer who can 
only access to the subregion A. Although we start off 
with a pure state with vanishing von Neumann entropy, 
the state corresponding to the reduced density matrix 
is generally a mixed state. p A contains the information 
on the quantum degrees of freedom traced out. The 
entanglement entropy is defined as the von Neumann 
entropy of the reduced density matrix, 



S A = -Trp A lnp A . 



(2.2) 
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If the system has a finite mass gap (or a finite correlation 
length), then the most of contributions to the entanglement 
between two regions comes from the field degrees of freedom 
near the boundary (see Fig. ||). That is, those in A and B sep- 
arated by more than the correlation length have no quantum 
correlations and do not contribute to entanglement. Accord- 
ingly, the entanglement entropy saturates above some critical 



length in quantum field theories with mass gap [10]. For in- 
stance, the entanglement entropy of the ground state of the 
Ising chain model in the non critical regime saturates at large 
/, Sa(1) — * c/3 log (%/a). Here % is the correlation length of 
the system and a the lattice spacing, c the central charge. 

Although an analytic proof is still lacking, numerical ev- 
idence of the existence of the mass gap in Yang-Mills theory 
have been accumulated in 30 years, and we expect that such 
a saturation can be seen in lattice QCD simulations. Some 
other properties of the entanglement entropy can be found in 

iin. 
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Figure 3: Schematic figure of 
spin system with a finite correlation 
length £, . Spin degrees of freedom 
in two regions A and B separated 
more than the correlation length do 
not have quantum correlations, and 
do not contribute to quantum entan- 
glement. 



3. Replica trick 



B 



In order to evaluate the entanglement entropy, we ap- 
plied the replica trick. The detail of the derivation is 
given in [O]. The point is that the entanglement en- 



tropy defined in Eq. ( |2.2[ ) can be represented in the form, 
Sa = — lim n _»i d /dn\rilrAPl- 

The trace of the rc-fh power of the reduced density 
matrix p A is given by the ratio of the partition functions, 



Tr p" A =Z(l,n)/Z n . 



(3.1) 



Here Z(l,n) is the partition function of the system having 
special topology, the «-sheeted Riemann surface. The field 
variables in the region A is periodically identified with the 
interval n/T (T is the lattice extent in the temporal direc- 
tion, corresponding to temperature) while in the region B 
the periodic boundary condition is imposed with the period l/T. The case for n = 2 is illustrated 
in Fig. ^. Z is the partition function with no cut. 
The entanglement entropy is then given by 

z(i, n y 



Figure 4: Schematic picture for the 
system with two cuts in x — t plane. In 
the region A (B), the periodic bound- 
ary condition is imposed with the pe- 
riod 2/T{X/T). 
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(3.2) 



The derivative of Sa (0 with respect to /, which is free of the ultraviolet divergence, can be expressed 
as follows; 



dS A (l) 



dl 



d_ 
dl 



lim In 

«->i on 



Z(l,n) 
Z n 



n—>i dl on 



(3-3) 
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That is, in order to calculate dS^/d/, we first evaluate the free energy of the system having n cuts 
with the length / of the cut, then take the derivative with respect to n and /, and take the limit 
n — » 1 . Thus, the evaluation of the entanglement entropy is reduced to calculate the free energy of 
the system with n cuts. 



4. Lattice setup and observables 



In numerical simulations, the derivative in Eq. ( p.3| ) have to be replaced by the finite difference, 
and we estimate the derivative by 



Km——F[A,n] 
«->! dl dn 



d , , , , ,. F\l + a,n = 21 — F\l,n = 21 

-lim F /,n + l ]-F[l,n }) ^ l - [ - l -. (4.1) 

dl n->l a 



In the first line, the derivative with respect to n is replaced by the finite difference between the free 
energies for n and n + 1 cuts. To go to the second line, we substitute n = 1 to the free energies. At 
this point, dF[l,n = l]/dl drops out since F[l,n = 1] does not depend on I. 

The differences of free energies can be evaluated numerically by introducing an 'interpo- 
lating action' which interpolate two actions corresponding to two free energies ^ |9|], Sj nt = 
(1 — a)Si[U] + aSi +a [U]. Si and Si +a represents the actions corresponding to F[l,n = 2] and 
F [I + a, n = 2} in Eq. ( |4.1[ ). It is easy to show that 

F[l + a,n = 2]-F[l,n = 2} = - da— lnZ(/,a) = / da (S, +a [(j)} - S,[(j)}) a . (4.2) 

Jo da Jo 

Here (■)« refers to the Monte Carlo average with the interpolating action (1 — a)Si[U] + aSi +a [U]. 
Therefore, the entanglement entropy can be evaluated numerically by updating gauge configura- 
tions with the interpolating action and calculating the action differences for various a and perform 
a numerical integration over a. In order to evaluate the integral in Eq. ([0|), we calculated the 
action differences from a = to 1 by the step 0.1, and employed the Simpson's rule to evaluate the 
integration numerically, which interpolates neighboring points by a quadratic curve. 

The lattice configurations are generated by the heat-bath Monte Carlo technique with the stan- 
dard Wilson plaquette action. In our simulations, the first 5000 sweeps are discarded for thermal- 
ization, and the measurement has been done every 100 sweeps. The number of configurations for 
each j8 and lattice size is around 3000 to 8000. 

In Fig. HI we plotted the action differences (5/=2 — 
5/_i ) on 16 3 x 32 at jS = 6.0. As is clear from the figure, 
the line connecting the data points crosses zero at about 
a = 0.5. Thus, most of the contribution cancels in inte- 
gration over a from to 1 though the absolute values of «- 
the differences are large at both end points. 



5. Simulation results 
Entanglement entropy 

The derivative of Sa(1) with respect to / is plotted 
in Fig. |6| dSA(l)/dl is normalized by the area of the 
common boundary, \dA\. 



Figure 5: The difference (5/ = 2 — S/=i) on 
16 3 x 32 at [i = 6.0. The integration from 
a = to a = 1 gives BSa/SI at I = 3a/2, 
the midpoint between I = a and / = 2a. 
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Figure 6: The derivative of the entanglement entropy Sa normalized by the area dA of the common 
boundary with respect to the length of the region A. The dashed curve is the fit of the data by the function 
c/l a with the fitted values c — 0.149(48), a = 3.06(20). The right panel shows the zoom up of the left panel 
to make near-zero region more visible. 

We observe that data on 12 3 x 24 and 16 3 x 32 agree within statistical errors. This implies 
that the derivative of the entanglement entropy is proportional to the area of the boundary as is 
expected. 

As is explained above, the entanglement entropy is closely related to the correlation length E, 
of the Hamiltonian. Thus, in small / regions, the entanglement entropy is expected to scale as I /I 2 
from the dimensional analysis. That is, BSa/SI behaves as l// 3 at small /. This behavior is exactly 
what the entanglement entropy in conformal field theory in (3+l)-dimensional spacetime shows. 
In order to confirm this, we fitted data with the function dSA/dl = c(l/l) a , and we obtain c = 
c = 0.149(48), a = 3.06(20), % 2 /ndf = 0.192. The fitted function is plotted in Fig. | by dasshed 
curve. 

The string feature of our result is that the derivative of the entanglement entropy hits zero at 
about I* = 0.6 ~ 0.7 [fm] and vanishes above. This means that the entanglement entropy does not 
increase with increasing I at large / and the QCD vacuum has a finite correlation length (or a finite 
mass gap). We note that the critical temperature T c of SU(3) pure Yang-Mills theory is estimated 
from the behavior of the Polyakov loop susceptibility as 280 [MeV], and i/T c ~ 0.714 [fm] fll3|]. 
This value and our result promise the identification, /* = l/T c . In other words, the critical length 
of entanglement entropy and (the inverse of) the critical temperature of the deconfinement phase 
transition may be identified. 
Entropic C-function 

The holographic analysis of the entanglement entropy for confining backgrounds have revealed 
that there is a transition from the connected solution to the disconnected solution corresponding to 
0{N 2 ) and 0(1) solutions, respectively [Q, |p. This predicts a jump in the entropic C-function, 
C(Z) = I 3 /\dA\dSA/dl, at some critical length I*, above which the C-function vanishes. This may 
be parallel to the confinement/deconfinement phase transition. 

The numerical result for the entropic C-function is shown in Fig. [7| We observe that the C- 
f unction takes non-zero value below I* ~ 0.6 [fm]. Above 0.6 [fm], the numerical data suffer from 
huge statistical errors. This is because BSa/SI is very small in this region while its statistical errors 
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Figure 7: The entropic C-function normalized by an area of the boundary, C(Z) = m-j^. 

do not so much depend on /. Thus, the relative statistical error becomes quite large at large /, and 
we cannot specify the critical length precisely. 

6. Summary and conclusion 

We studied the entanglement entropy of the QCD vacuum in SU(3) Yang-Mills theory using 
lattice Monte Carlo simulations. The entanglement entropy is defined as the von Neumann entropy 
of the reduced density matrix which is obtained by tracing out the degrees of freedom in one of two 
complementary regions. It measures the quantum correlation between subregions. We find that the 
entanglement entropy scales as 1/Z 2 at small / as in the conformal field theory. The derivative of 
the entanglement entropy with respect to Z hits zero at about I* = 0.6 ~ 0.7 [fm] and it vanishes 
above this length. It implies that the Yang-Mills theory has the mass gap of the order of 1/Z*. This 
value is very close to the critical temperature of the confinement/deconfinement phase transition. 
The entropic C-function is suffered from large statistical errors and we cannot specify the critical 
length precisely. Simulations with the renormalization-group improved action will improve this 
situation. 
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A. Riemann sheets and gauge invariance 

In lattice gauge theory, the Monte Carlo update of link variables on rc-sheeted Riemann surface 
needs special care in order not to violate the gauge invariance. 

Consider a lattice of the volume 8 3 x 4 and double the lattice in the temporal direction to 
calculate the free energy of the lattice with two cuts, see Fig. |8|. We denote links on two lattices 
as {U } and {U'}. The periodic boundary condition in the region B is such that (for simplicity, we 
suppress the Lorentz indices) 

• U (it, / = 0) = U(x,t = 4) (red arrows in Fig. g) 

• U'(x,t = 0) = U'(x,t = 4) (blue arrows in Fig. |), 
and that in the region A is 
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Figure 8: (Left) Riemann surface structure of the lattice system with two cuts. Links with same colors 
indicate periodically identified links. (Right) Two Riemann sheets projected on a single plane. In order to 
update the link variable denoted by U without violating the gauge invariance, we have to use the products of 
links (plaquettes) enclosing the area a and c, not a and b. 

• U(x,t = 0) = U'(x,t = 4) (green arrows in Fig. ||) 

• U(x,t =4) = U'(x,t = 0) (magenta arrows in Fig. [8]). 

In order to update the link variable denoted by U in Fig. |8j, we need six plaquettes, two of them are 
lying inx — t plane. Naively, we use the plaquettes denoted by a and b in the right panel for two 
of six plaquettes to update the link U. However, this leads to the violation of the gauge invariance 
because the plaquettes a and b belongs to different Riemann sheets. Therefore, the link U lies in 
the upper Riemann sheet in Fig. || does not make a closed loop with the staple b, meaning the loss 
of gauge invariance. We have to use the plaquette denoted by c to update the link U otherwise the 
gauge invariance is not preserved. 

References 

[1] L. Amico, R. Fazio, A. Osterloh and V. Vedral, Rev. Mod. Phys. 80, 517 (2008), [quant-ph/0703044]. 

[2] S. Ryu and T. Takayanagi, Phys.Rev.Lett. 96, 181602 (2006), [hep-th/0603001]. 

[3] T. Nishioka, S. Ryu and T. Takayanagi, 0905.0932. 

[4] I. R. Klebanov, D. Kutasov and A. Murugan, Nucl.Phys.B 796, 274 (2008), [0709.2140]. 

[5] T. Nishioka and T. Takayanagi, JHEP 0701, 090 (2007), [hep-th/061 1035]. 

[6] A. Velytsky, Phys. Rev. D77, 085021 (2008), [0801.4111]. 

[7] P. V. Buividovich and M. I. Polikarpov, Nucl. Phys. B802, 458 (2008), [0802.4247]. 

[8] G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, 0710.4197vl. 

[9] Z. Fodor, PoS LAT2007, 011 (2007), [0711.0336]. 

[10] M. Wolf, F. Verstraete, M. Hastings and J. Cirac, Phys.Rev.Lett. 100, 070502 (2008), [0704.3906v2]. 

[11] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge 
University Press, 2000). 

[12] P. Calabrese and J. L. Cardy, J. Stat. Mech. 0406, P002 (2004), [hep-th/0405152]. 

[13] Y. Iwasaki, K. Kanaya, T. Kaneko and T. Yoshie, Phys. Rev. D56, 151 (1997), [hep-lat/9610023]. 



8 



